function lower_bound = gmmBayesLowerboundVerbeek(mix, data, vars, reg)


% compute sufficient statistics
[mix.sample_means mix.sample_covars] = computeSufficientStatisticsVerbeek(data, vars, mix.covar_type, reg);


% E[ln p(Y|C, X=Y)]
q1 = 0;     % not actually zero, just never changes as the posterior is observed

% E[ln p(X|Z, mu, Lamba)]
q2 = ElnpX_ZmuLambda(mix, vars);

% E[ln p(Z|pi)]
q3 = ElnpZ_pi(mix, vars);

% E[ln p(pi)]
q4 = Elnppi(mix);

% E[ln p(mu, Lambda)]
q5 = ElnpmuLambda(mix);

% E[ln q(X=Y|Z)]
q6 = 0;     % not actually zero, just never changes as the posterior is observed

% E[ln q(Z)]
q7 = ElnqZ(vars);

% E[ln q(pi)]
q8 = Elnqpi(mix);

% E[ln q(mu, Lambda)]
q9 = ElnqmuLambda(mix);


lower_bound = q1 + q2 + q3 + q4 + q5 - q6 - q7 - q8 - q9;

% fprintf('%.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g %.10g\n\n', [q1 q2 q3 q4 q5 -q6 -q7 -q8 -q9]);
